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Summary: Semi-parametric methods are often used for the estimation of intervention effects on 
correlated outcomes in cluster-randomized trials (CRTs). When outcomes are missing at random 
(MAR), Inverse Probability Weighted (IPW) methods incorporating baseline covariates can be used 
to deal with informative missingness. Also, augmented generalized estimating equations (AUG) 
correct for imbalance in baseline covariates but need to be extended for MAR outcomes. However, 
in the presence of interactions between treatment and baseline covariates, neither method alone 
produces consistent estimates for the marginal treatment effect if the model for interaction is not 
correctly specified. We propose an AUG-IPW estimator that weights by the inverse of the probability 
of being a complete case and allows different outcome models in each intervention arm. This estimator 
is doubly robust (DR), it gives correct estimates whether the missing data process or the outcome 
model is correctly specified. We consider the problem of covariate interference which arises when the 
outcome of an individual may depend on covariates of other individuals. When interfering covariates 
are not modeled, the DR property prevents bias as long as covariate interference is not present 
simultaneously for the outcome and the missingness. An R package is developed implementing the 
proposed method. An extensive simulation study and an application to a CRT of HIV risk reduction- 
intervention in South Africa illustrate the method. 

Key WORDS: Augmentation; Cluster-randomized trials; Generalized estimating equation (GEE); 
Interactions; Interference; Inverse probability weighting (IPW); Missing at random (MAR); Out¬ 
come Model; Propensity Score; R package; Semi-parametric methods. 
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1. Introduction 

In clustered randomized clinical trials (CRTs), the unit of treatment assignment is a cluster 
of subjects, which we also refer to as community. In such settings, outcomes are likely to 
be correlated among subjects within the same cluster. Often used for estimation, gener¬ 
alized estimating equations (GEE) based on semi-parametric methods (Zeger and Liang, 
1986) target marginal effects of treatment. Within clusters, dependence can be modeled 
using a working correlation structure. Compared to mixed effects models, this approach has 
the advantage of focusing on population average effects rather than cluster specific effects 
(which are equal for continuous outcomes) and requires fewer parametric assumptions on 
the outcome distribution (Hubbard et ah, 2010). Moreover, because both the outcome and 
the missing data mechanism can be modeled, this approach allows doubly robust estimation, 
which is impossible with mixed effect models. Finally, this approach to estimation is robust 
to misspecihcation of the correlation structure. However, challenges arise in developing a 
consistent and efficient estimator of marginal treatment effects; these include the need to 
adjust for missing data and accommodate covariate interference (wherein a subject’s outcome 
may be affected by covariates of other subjects) and interactions (wherein the effect of 
treatment varies by covariate-defined subgroups). We propose a method that addresses these 
issues and is practical to implement for evaluating novel interventions in CRTs. 

In CRTs, covariates may be fully observed even if the outcome is missing. When data are 
assumed missing completely at random (MCAR) - i.e. the observed process is independent 
of observed and unobserved information (Rubin, 1976) - the standard GEE approach pro¬ 
vides consistent and asymptotically normal (CAN) estimators. If the pattern of missingness 
depends on observed information but not on missing data, the data are said to be Missing 
at Random (MAR). In this setting, the standard GEE may yield biased estimates although 
likelihood-based approaches, such as mixed effect models, can provide unbiased estimators. 
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Imputation (Paik, 1997) or reweighing (Robins et ah, 1995) methods can correct for this bias. 

Although useful if the missingness mechanism is not completely known, multiple imputation 
requires correct specification of the joint distribution of the outcomes, which is especially 
difficult when they are correlated and the cluster sizes are large (Beunckens et ah, 2008). 

In this article, we consider the Inverse Probability Weighting (IPW) approach to analyze 
incomplete data. If the model for the missingness mechanism represents the MAR data 
generating process, the IPW estimation provides CAN estimators of treatment effects by 
reweighing complete cases according to the probability of being observed (Liang and Zeger, 

1986; Robins et ah, 1994). 

Recent methodological developments improve estimation efficiency by leveraging baseline 
covariates; they may be based on targeted maximum likelihood (Moore and van der Laan, 

2009) and on augmentation (Robins et ah, 1994; Robins, 2000; Tsiatis et ah, 2008; Zhang et ah, 

2008). Stephens et ah (2012) developed the augmented GEE (AUG) methods in the setting 
of dependent outcomes such as in CRTs. The AUG adds a term to the standard GEE 
which relates the outcome to covariates and treatment. Randomization assures that the 
AUG is CAN even in the case of OM misspecihcation. However in the case of outcome data 
that are MAR but not MCAR, the AUG may be biased. There exists theory for extending 
these methods to MAR data for individual randomized Trials (RTs) with possibly correlated 
data (Van der Laan and Robins, 2003; Glynn and Quinn, 2010), we focus on the details of 
implementing the methods in CRTs. 

The term interference can refer to different types of relationships among exposures, out¬ 
comes and covariates. Interference in RTs arises when one subject’s treatment may impact the 
outcomes of other subjects (Rosenbaum, 2007; Vansteelandt, 2007; Tchetgen Tchetgen and VanderWeele, 
2012; Hudgens and Halloran, 2012). A similar phenomenon, confounding by clusters, has 
been discussed in the context of observational studies (Seaman et ah, 2014); we will refer to 
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such confounding as exposure interference. In CRTs all subjects within a cluster receive the 
same treatment; hence if the clusters are independent as typically assumed in practice, there 
is no exposure interference measured at the cluster level. Therefore, any choice of working 
correlation structure for the standard GEE will give a consistent estimator of the marginal 
treatment effect (Pepe and Anderson, 1994). We will investigate covariate interference among 
individuals nested within clusters: the setting in which one subject’s covariate may impact 
the outcomes of other subjects. 

The IPW and the AUG can be combined in a doubly-robust method we refer to as the 
DR; we investigate its properties regarding robustness to misspecihcation of the missing data 
and outcome generating process. By considering a variety of data generating mechanisms, we 
investigate settings in which the DR has advantageous properties (consistency and precision) 
compared to the IPW and the AUG, and discuss the impact of covariate interference and 
treatment-covariate interactions. This paper is organized as follows. Section 2 introduces 
notation and assumptions for the IPW and the AUG GEE approaches. Section 3 describes the 
DR approach, investigates CAN properties and discusses the issue of covariate interference. 
Section 4 provides a motivating example with data arising from a CRT of an HIV / Sexually 
Transmitted Infection (STI) risk reduction intervention in South Africa (Jemmott III et ah, 
2014). Simulation studies regarding bias, relative efficiency and coverage are described in 
Section 5, and concluding remarks are made in Section 6. 


2. Notation, basic models and assumptions 

2.1 Notation for CRTs and marginal treatment effect 

We consider a study design in which a vector of P baseline covariates Xij = (Ah,... ,X?) 
and outcome Yjj are recorded for each subject j = 1,..., n* in community i = 1,..., M. 
The sample size within each community is assumed fixed by design and non-informative. 


4 


Biometrics, December 2015 


Our setting compares two arms (treated Aj — 1 and control At = 0); the probability of 
treatment assignment is known and given by p = P(Ai = 1); extension to a greater number 
of treatments is straightforward but complicates the notation. In this article, the outcome 
Yi = \Yij]j=i,...,ni is assumed to be continuous, but extension to other types of outcomes is 
straightforward. The vector R { = [Rij]j=i,..., ni is the indicator of missingness; Y tJ is observed 
when Rij = 1. The matrix of covariates X t = [Xij\j = i^.^ ni is assumed to be fully observed 
and consists only of pre-exposure covariates measured at baseline. 

Interest lies in estimating the marginal effect of the treatment given by = E(E(Y tJ \Ai = 
l,Xi) — E(Yij\Ai = 0,Xj)). For estimating M£, we make inference about the parameters 
(3 = (/3 0 ,(3 a ) t indexing the marginal model #(a%(/3, A)) = g( E { Y ij I A)) = Po + PaA, 
where /x, : (/3, Ai) = [/%(/3, A)]j=and g is a one-to-one link function, which is an identity 
function in this article. Of particular interest, Pa is equal to M^. Of note, extension to binary 
outcome Y t j using a logistic function for g and considering odd-ratios is based on the same 
reasoning. 

When the outcome is believed to be MCAR, the missingness process is independent of Xj, 
Ai , and Y^. If one assumes MAR and the missingness pattern is monotone, the probability of 
missingness can be estimated by a multistep approach by decomposing a monotone missing 
pattern into multiple uniform missing data models (Robins et ah, 1994; Li et ah, 2011). In 
CRTs, any component of Y\ can be missing; hence the missingness pattern is non-monotone. 
Therefore, we make a stronger assumption than MAR that we refer to as restricted MAR 
(rMAR): the probability that the outcome for one individual is missing is independent of 
all outcomes in the cluster, conditional on baseline exposure Ai and cluster characteristics 
X[. The conditional probability that the outcome is observed is denoted 7Tij(Xi,Ai) = 
P(Rij = 1| Xi,Ai) and is called the propensity score (PS). When data are rMAR, ignoring 
missing data leads to biased inference if missingness depends both on X; and Ai. This is 
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because the presence of missing data no longer assures balance on average of confounding 
factors between treatment arms. Therefore, analysis must include adjustment for missing 
data; appropriate models for this adjustment may require treatment-covariate interactions, 
which may be difficult to specify and require many parameters. Combining the IPW and 
the AUG, which this paper proposes, makes it possible to obtain consistent estimates of 
the marginal effect of treatment without explicitly specifying interaction terms while also 
improving efficiency. 

2.2 Inverse Probability Weighted Generalized Estimating Equations (IPW) 

In order to account for missing data, semi-parametric estimators based on the IPW are found 
by solving the estimating equation 1: 

M 

0 = MA,)}, ( 1 ) 

---' 

^(XuRiAuP, v w ) 

where Di = 1 is the design matrix, Vi is the covariance matrix equal to U 1 / 2 C(ct)U ]/ 2 

with Ui a diagonal matrix with elements var(r/jj) and C(ot) is the working correlation 
structure with non-diagonal terms ot. For example, for an independence correlation structure 
ct are zero; for exchangeable all the elements of ot are identical. Parameters ct could also 
depend on the treatment group C(a(A*)) but we do not consider this possibility in our 
implementation. In this article, we estimate the ct parameters using moment estimators 
from the Pearson residuals as in McDaniel et al. (2013). The n i x rq matrix of weights 
is Wi(Xi, A i: r) w ) = diag[Ri j /'n- ij (X i ,A il ri w )\ j=1 ^ ni , where the PS is derived by fitting 

a binary response model to the indicator Rij regressed on Ai and a subset of X, - say 

using a logistic regression. The rj w are nuisance parameters estimated in the PS. A nec¬ 
essary assumption for this method is that probabilities for the PS are bounded away from 
zero. Several authors have noted the instability that may arise from small probabilities of 
observation (i.e. large weights) and proposed use of stabilized or truncated weights; see 
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Seaman and White (2013) for a review. To ensure that the IPW provides a CAN estimator, 
the PS must include all covariates that are associated simultaneously with the missingness 
and outcome (Brookhart et ah, 2006), including those that involve interaction with treatment 
(Belitser et ah, 2011). 

2.3 Augmented Generalized Estimating Equations (AUG) 

For settings with complete data, Stephens et ah (2012) proposed the AUG estimator which 
can improve efficiency relative to the standard GEE by incorporating baseline covariates. 
The AUG is constructed by subtracting from the set of GEEs the orthogonal projection 
of the standard estimating function onto the span of scores corresponding to all smooth 
parametric models for the treatment assignment mechanism given covariates. The AUG is 
given in Equation 2: 



0 


il>i(Yi,Aa 3 ) 


+ Y,P a ( 1 -P) 1 ~ aD T v 7 1 (B(X i ,A i = a,T tB )-fji i (0,A i = aj) . ( 2 ) 


a=0,l 


The term YAi, (3) is similar to iPi(Yi, Ri , Ai, /3, r] w ) in Equation 1 for the IPW except 
that Wi is set to identity because there is no adjustment for missing data. Definitions for 
Dj and Vi remain the same. The function B(Xi,Ai = a,r) B ) is an arbitrary function of 
X, given for each treatment arm. The r] B are nuisance parameters that must be estimated. 
The estimator in Equation 2 is most efficient if B(Xi, A t = a , rj B ) models the outcome and 
is equal to E(Yij\Xi,Ai = a) (Robins et ah, 1994; Zhang et ah, 2008). Thus, B(X i) A i = 
a, rj B ) is called the outcome model (OM). In the absence of missing data, the AUG remains 
consistent even if the OM is not correctly specified (B(X i} Ai = a,rj B ) ^ E(Y ij \X i ,A i = 
a)). Correct specification can lead to substantial efficiency gains compared to the standard 
GEE. Moreover, in presence of treatment-covariate interactions, it is useful to fit a different 
regression model for the OM for each treatment group, e.g. B(Xi,Ai = a,rj B ) = 7o + 
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7 r^-lj with rj B = ( 7 °,..., 7 p, , 7 p), thereby obviating the need to fit covariate- 

treatment interactions terms. In presence of rMAR, the AUG does not ensure consistent 
estimation; instead, one must combine the AUG with the IPW as we show below. 


3. Methods to accommodate missing data, treatment-covariate interactions 
and covariate interference in CRTs 

3.1 Doubly Robust Augmented IPW Generalized Estimating Equations (DR) 

We extend the AUG in Equation 2 to account for missing data using the IPW in Equation 
1 by subtracting from the set of GEEs the orthogonal projection of , ij) i (Y i , Ri, (3,rj w ) 
onto the span of scores corresponding to all smooth parametric models for the missing data 
process and the treatment assignment mechanism given covariates (Tsiatis, 2006). This gives 

the following estimating equation (see Web-Supplementary Material B for details): 

m r 


0 


E 


DjV^W^Xi, At, lw) (Yi 


B(Xi,Ai,rj B )) 


+ E P a{l ~ Py~ aD i V i 1 Ai = a ’ Vb) - PiO 3, A = a )) , 

a=0,l 


(3) 


M 

= i, Ri, Ai, Xi, (3, rj W) rf B ). 

2—1 

The Di, Vi and the PS are defined such as in Equation 1, the OM denoted B(Xi,Ai = 
a, rj B ) is defined for each treatment group such as in Equation 2. The estimator denoted 
Aug is found by solving the estimating equation given in equation 3. Although analytic 
solutions sometimes exist, coefficient estimates are generally obtained using an iterative 
procedure such as the Newton-Raphson method. To get f3 aug we use the estimated PS 
(Wi(Xi,Ai,rj w )) and estimated OM (. B(Xi , Ai,rj B ))- A mentioned above, treatment- 
covariate interactions can be accounted for by fitting OM regressions separately by treatment 
group. One could also estimate parameters of the PS model separately by treatment groups. 
This approach, however, may provide less stable results due to variability in the calculation 
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of weights. Iii this paper, fj w in Wj(Xj, A i: rj w ) are obtained using a logistic regression and 
f) B in B(Xi, Ai,f] B ) are obtained using a linear regression. Thus, we treat R t] and R tj > as 
conditionally independent given Ai and X t . In the presence of correlation of Rij and Riji , 
one might be able to improve efficiency of estimation of 7pj and therefore of the marginal 
treatment effect by accounting for this correlation. Of note, estimation procedures other 
than generalized linear models could also be used to compute the OM and the PS values. 
The DR estimator is doubly robust in the sense that it is CAN under correct specification 
of either the OM (i.e. B(X i: A^tj B ) = EfYifAi, Xf)) or the PS (i.e. iTij(Xi, A^f} w ) = 
P{Rij = l\Xi,Ai)) (see Web-Supplementary Material Section Cl). Implementation in R is 
available on the CRAN in the package ’CRTgeeDR’. Source code had been made available 
as Web-Supplementary material. We note that in contrast with several existing software 
packages (for example proc GENMOD in SAS (2015)), our implementation of the weighted 
GEE, which uses V^W^X,, A h i lw ) instead of W] , 2 (X h A t , ■n w )V^ 1 W 1 / 2 (X i , A h rj w ), 


guarantees consistency for all choices of working correlation structure (see details in Web- 
Supplementary Material Section C2 and D). 

3.2 Variance of the DR estimator 

The variance of (3 aug is estimated by the sandwich variance estimator. There are two external 
sources of variability that need to be accounted for: estimation of rj w for the PS and of r] B 
for the OM. We denote VI = ((3, rj w , rj B ) the estimated parameters of interest and nuisance 
parameters. We can stack estimating functions and score functions for VI: 




1 &i(Yi,Xi,Ai,(3,r) w ,ri B ) ^ 
SV(Xi,Ai, r] w ) 

y Si {Xi, Ai,rj B ) J 


where SY and Sf represent the score equations for patients in cluster i for the estimation 
of rj w and r] B in the PS and the OM. A standard Taylor expansion paired with Slutzky’s 
theorem and the central limit theorem provide the sandwich estimator adjusted for nuisance 
parameters estimation in the OM and PS. We refer to this as the nuisance-adjusted sandwich 
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estimator: 


Var{Et) = E 


dUi( O) 

on 


i -i T 


E [Ui(n)Uj(Q)\ E 


dUi{Q) 

on 


-i 


(4) 


*adj 


n - 1 

adj 


The variance estimator var((3 ) is obtained by estimating unknown quantities upon sub¬ 


stituting empirical means for expectations and ft = for f2. Thus, the term A ad j 

is given by A Yf=i U i{H)U i{VL) T and T adj is given by ^ Y^U 
In small sample settings, it is likely that this estimator of the variance of 0 aug is bi¬ 


ased. We implemented Fay’s bias-correction approach, which is particularly suitable for 


M-estimators (Fay et al. 2001). The term A ad j in Equation 4 is replaced by A f ay given 


l i \—\M 

b y m Ei=i 


H; 


\ T ' 


U t (n) iHiUiiQ)) 


H 


i [jj] 


1 — min(q, 


' dUi(Sl) r i \ 

' dfl 1 adj) [jj] 


, where H, is a diagonal matrix with diagonal terms 
, q = 0.75 is a frequently-used bound. 


3.3 Definition of covariate interference and implication for analysis 

In previous sections, we discussed covariates measured on the index subject (j), but other 
subjects’ (j r ) covariates may also impact the outcome for the index subject. An example of 
a potentially interfering covariate is described by Kaiser et al. (2011) who found a positive 
association between age of partner and infection with HIV. Similarly, the characteristics 
of subgroups to which the index case belongs (household, neighborhoods, ...), whether 
known or not, may be interfering covariates (Brumback and He, 2011). In this paper, we 
consider the phenomenon of covariate interference where there exists at least one individual 
j' j such that E(Yij\Xij) ^ E(Yj.j\Xij, Xij/), where X,j represent the vector of all 
measured baseline covariates. That is, even after all covariates for the index subject j have 
been included in the model, the covariates of individuals other than the index subject 
still affect the outcome of the index subject j; we refer to such covariates as interfer¬ 
ing covariates. See Pepe and Anderson (1994) for a similar definition in longitudinal data 
and see Seaman et al. (2014); Liu and Hudgens (2014) for an analogous definition in non- 
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randomized clustered data in the context of confounding by cluster and interference. Refer to 
Web-Supplementary Material Section A for a causal interpretation of covariate-interference. 

When interfering covariates affect either the outcome (E( y Y i j\X i j) ^ E ( Y VJ \ X l3 , X,ji )) or 
the missingness process (E(R i j\X i j) E{Rij\X i3 , X i3 >)), but not both, the DR estimator is 
CAN even if the interfering covariates are not included in the models, provided that either 
the PS (P(Rjj\Xij, Ai)) or the OM (E(Yij\X^, Ai = a)) is correctly specified; that is, either 
the PS or the OM includes all the covariates X %3 involved in the same functional form as 
in the data generation processes. Accounting for covariates interference in the OM increases 
efficiency if and only if they predict the outcome. When interfering covariates impact both 
the outcome and the missing data generating processes, they must be included in either 
the OM or the PS models in a way that correctly represents the data generation processes. 
Thus, it will ensure that the DR estimator will be CAN if a correct model for either the OM 
(E(Yij\Xi, Ai = a)) or the PS (P(R i j\Xi,A i )) is specified, where the X i3 are replaced with 
Xi in the formulas above. We acknowledge that this model for interfering covariates is not 
likely to be known and can be difficult to identify. Different cluster sizes and sub-clustering 
structures (such as households) may make infeasible the use of regression techniques in the 
OM or the PS because of the potentially different dimensions of the individual and interfering 
covariates. Cluster summary measures such as the mean or maximum of individual covariates 
in the cluster (or sub-groups in each cluster) may nonetheless be useful in incorporating 
interference covariates in models (Brumback et ah, 2010). 


4. Application 

4.1 Description of the SAM study 

We analyze data from the “South African Men” (SAM) study which randomized 22 pair- 
matched clusters to a health-promotion intervention (control) and an HIV /STI risk-reduction 
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intervention in a CRT design; the study included 1181 South African men who have sex with 
women. A complete description of the study design can be found in (Jemmott III et ah, 
2014). We focus on a cross-sectional analysis of these data after one year and ignore matching. 
The primary outcome of our analysis is the overall percentage of acts of protected intercourse 
among the total number of acts of intercourse. When the total number of acts of intercourse 
is zero, we set the percentage to 100%, as no exposure implies no risk. Secondary outcomes 
are the percentages of protected acts of intercourse by type of partnership and type of 
intercourse (vaginal and anal sex with main and casual partners). Descriptive statistics 
for these outcomes, including proportion of missing observations by type of partner and 
intercourse are provided in Table 1. Slightly more observations are missing in the HIV/STI 
intervention group (20.8% versus 17.5%). The overall protection percentage after one year 
are about 64% for the HIV/STI intervention compared to 60% for the control group. 

As the proportion of missing baseline covariates was less than 0.1%, we consider them to be 
MCAR and exclude observation with missing covariates from the analysis. No community 
sub-structure, such as household or neighborhood structures, was described in the SAM 
study. Here we consider potential interfering covariates at a cluster level by taking the 
mean (or mode for qualitative variables) of baseline measures in the community: X/ = 
/- X/ ? =i n , ■ F° r example Hawkes et al. (2013) demonstrated that the mean religiosity 
score for a community, defined as the mean of individual religiosity score in the community, 
may have an impact on each individual outcome and missingness in particular regarding 
sexual behaviors. Table 1 describes socio-demographical individual-level variables and inter¬ 
fering covariates. We provide p-values for Wald tests testing the association of covariates 
and treatment-covariate interactions with the outcome and the missingness indicator. In 
this study, there is evidence of interactions of individual covariates with treatment for both 
the outcome and the missing data generation processes. However, the interfering covariates 
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defined here do not appear to be significantly associated with both the outcome and the 
missing data generation process. 

[Table 1 about here.] 


4.2 Results 

We analyze these data with the GEE, the AUG, the IPW and the DR using both indepen¬ 
dence (-1) and exchangeable (-E) working correlation structures. Variables for the PS, and the 
OM were selected using a forward stepwise regression (separately for each treatment group) 
from among all the individual covariates X,j presented in Table 1. We did not include the 
interfering covariates (X ? .) in the analysis as none impacted both outcome and missingness 
processes (Table 1). We used the step function in R based on the AIC criterion. Results of 
these selections are given in Web-Supplementary Material F. We describe here the results for 
the primary outcome. The amount of missingness is larger in the treated arm and increases 
with age; it decreases with religiosity, good health score, and exercise. The OM patterns are 
substantially different for treated and control; the only common variable is the CAGE score. 
In both arms lower alcohol consumption is associated with a greater percentage of protected 
acts of intercourse. Results are presented in Table 2 for primary and secondary outcomes. 
With the DR-E, we observe a significant difference of 7.4% (sd=2.9%, p=0.01) in the overall 
percentage of protected intercourse in the ffIV/STI intervention group compared to the 
control group. Analyses of the secondary outcomes suggest that this result is mainly driven 
by condom use during vaginal intercourse with a marital partner. The HIV/ST1 intervention 
has no significant impact on other outcomes. Using the DR rather than the standard GEE 
or the AUG has an impact on the treatment effect estimates and associated standard errors 
(SE). The difference between these approaches is apparent in the magnitude and direction 
of the marginal treatment effect estimate. For example, the analysis for the GEE-1 (3.8 [- 
1.0; 8.5]) does not demonstrate a significant effect of the fflV/STl intervention on overall 
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percentage of protected intercourse, whereas this effect is stronger and significant for the 
DR-1 (7.3 [1.6; 13.0]). Both the GEE-I and the AUG-I (5.4 [2.2; 8.7]) are probably biased 
due to missing data. Using the DR instead of the IPW leads to an increased magnitude of 
the treatment effect and an increased level of statistical significance: for example, the DR-E 
(7.4 [1.73; 13.0]) compared to the IPW-E (3.4 [-1.4; 8.3]). 


[Table 2 about here.] 


5. Simulation Studies 

5.1 Properties of the DR estimator 

We consider a setting with continuous outcome Y and assignment of treatment Ai at a 
cluster level with probability p — 1/2. We generate a normally distributed covariate XIij 
(independent of Ai) with mean 1 and a standard deviation of 5. For each individual, we 
define a covariate XI j. which is the mean of XI for all the subjects in the same cluster: 



X3 i, are defined as was XI j. and are possible interfering covariates. The model for simulation 
is given in Equation 5: 



The parameters f3° = (/3(/, /3®, /3p, /3^) are the regressors associated with intercept, 


treatment, covariate, interfering covariate, treatment-covariate interaction for the outcome 


model. Parameters (3 M are the same for the missing data generating process. Scenarios with 


low correlation among cluster (0.05) were simulated with ef ~ A/”(0, 0.05) and ~ W(0,1.0) 
for cluster and individual random errors; scenarios with high correlation (0.2) were simulated 
with ef ~ A^(0,0.25) and ~ AT(0,1.0). True correlation structure is exchangeable. We 
investigate small sample (M = 10 and rq = (10, 20, 30) with probability 1/3 each) and large 
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sample (M = 100 and rii = (90,100,110) with probability 1/3 each) properties. In each 
scenario, we generate 1000 replicates of datasets. 

We evaluate the double robustness of the DR estimator in the setting of large and small 
sample with low correlation, but similar results are observed for large correlation. We inves¬ 
tigate models of analysis with OM and PS correctly specified (TRUE), misspecihed (MISS) 
and partially specified omitting treatment-covariate interactions (NONE). Table 3describes 
the data generation process, provides the formulations of the models of analysis, and shows 
the results from analysis; on average, 26% of outcomes were missing and the average ICC 
was 0.08. When there is no missing data, the traditional GEE is consistent because of 
randomization. When outcome data are MAR but not MCAR, the GEE and the AUG 
analysis are biased (-1.7 for the GEE-I and -1.8 for the AUG-I). When either the OM or 
the PS models or both are correctly specihed there is negligible estimated bias for the DR 
- a finding that confirm consistency. In small samples, this bias is bigger when only the PS 
is correct because the weights are estimated with lower accuracy. Using the more common 
choice of implementation for the weighted GEE W l J 2 {xi w )V^ l W l J 2 (Vw) leads to very high 
bias if an exchangeable correlation structure is used (0.374 if the OM is correct and 858 if it is 
not, for large sample). When the OM is correct the coverage remains around 95% (see Table 
2 in Web-Supplementary Material E). Using V’^ 1 Wi(X i , A,, rj w ) in the implementation 
of weights addresses this problem and permits the use of correlation structures other than 
independence. The IPW with correct PS also corrects the bias (-0.01) but is less efficient 
than the DR approach; coverage is close to the nominal value of 95%. In small samples, 
the empirical SE are underestimated. By contrast, in the large sample setting, using the 
nuisance-adjusted sandwich estimator for the DR leads to good estimates of the asymptotic 
SE (0.0263) compared to the empirical SE (0.0266) over 1000 replicates. Moreover, we observe 
that the coverage using the DR is comparable to that of the GEE with complete data. 
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Finally, we note that when the treatment-covariate interactions are ignored in the PS and 
only accounted for in the OM by fitting a different regression in each treatment group, the 
DR approach is also consistent and achieve same precision as when both the PS and the 
OM are correct (0.0014 and SE=0.027 for OM.TRUE.PS.NONE and 0.0013 SE=0.029 for 
OM.TRUE.PS.TRUE). 

[Table 3 about here.] 

Table 4 presents the results of analyses with the GEE, the IPW, the AUG and the DR 
that investigate the impact of correlation of the outcome in the data with small and large 
sample. The average percentage of missing outcomes is 23%; the average ICC is 0.04 for 
low correlation and 0.21 for high correlation. We analyzed the data using a PS and an OM 
model that was fit using a stepwise variable selection from among all of the individual and 
interfering covariates described above. The GEE and the AUG estimates are systematically 
biased because there is no correction for missing data. The IPW is also biased because the PS 
is incorrect in that it omits treatment-covariate interactions. The DR estimates are consistent 
in all analyses. In small sample settings, the empirical SE is underestimated even when 
using nuisance-adjusted SE, but estimation is improved by Fay’s correction. Nonetheless, 
the coverage remained lower than 86%, but it improves for large samples. Finally, when 
there is low correlation in the outcome, the robust SE better approximate the empirical SE. 

[Table 4 about here.] 

5.2 Simulations mimicking the SAM Study 

To consider more complex settings, we mimic the SAM study (see Section 4). We simulate 
the following individual-level covariates: employment (EMP ~ 13(0.25)), marital status 
(MARI ~ 73(0.23)), age (AGE ~ J\f( 27; 7)), religiosity (REL ~ W(0, 0.8)), the CAGE score 
(from a multinomial of probabilities CAGE ~ A4(0.3; 0.1; 0.1; 0.2; 0.3) for modalities 0,1,2,3 
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and 4), the HIV score (HIV ~ A/"(14; 4)) and the condom knowledge score (CDM ~ A^(3; 1)). 
Interfering covariates are generated as means for quantitative variables or modes for qualita¬ 
tive variables of the individual-level variables in each of the community (as was done for XI, 
X2 and X3 in Section 5.1). We generate data from the model in Equation 6. In simulating 
the outcome, we add cluster random errors to create an exchangeable correlation structure 
with ef ~ A/”(0, 5) and an individual random effects ~ A/"(0,4). This provides an outcome 
correlation among clusters of 0.07. We analyzed the data using a PS and an OM composed of 
all the covariates described above with a stepwise variable selection. Table 5 shows the bias, 
SE, and coverage of the methods we consider based on 1000 replicates for the estimation 
of the parameter M^ = 5.73. The percentage of missing outcomes is 21% and the average 
empirical ICC is 0.06. 

Yij = 60+40Ai—9.0EMP;j—S.OMARIij+l.OCDMy+S.ORELij 

T Ai[—2.0AGEij+8.5EMPij+3.5MARIij+1.5HIVij—2.0CAGEi 3 '+2.0RELi 3 ] 

v -v-' 

Interactions 

r 0.5AGE~-7.0CDMi.-5REL~+1.0HrvT +£f+eg. ^ 

covariate interference 

logit[P(Rij= 0)] = —3.0+2.0Aj+0.01AGEij—0.1HIVi 3 '+Ai[—0.1AGEjj—0.2HIVjj] 

V -V-' 

Interactions 

+ 0.02AGE7+0.2CDM7+0.2CAGE7 

V ---V---' 

covariate interference 

Table 5 provides the estimates the marginal treatment effect for small sample and for the 
same sample size as that of the SAM data. The GEE, the AUG and the IPW yield biased 
results whereas the DR has small bias justifying its use to analyse the data even ignoring 
covariate interference. Fay’s correction with coverage around 92% in small sample and 95% in 
large sample achieve good accuracy. Figure 2 in Web-Supplementary Material C3 represents 
the histograms of estimates over the 1000 replicates together with the true value of marginal 
treatment effect. It displays the bias of the GEE, the AUG and the IPW estimators compared 
to the DR and supports the approximate normal distribution of the DR estimator. 
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[Table 5 about here.] 


6. Discussion 

We propose a doubly robust method for the estimation of the marginal effect of treatment in 
CRTs with continuous data subject to rMAR - an assumption that arises because missingness 
is non-monotone in CRTs. Extension to binary or other outcomes is straightforward, provided 
that there is a one-to-one link function h such that: ffij = h(Xi,Ai). We extend the IPW 
approach proposed by Robins et al. (1995) and the AUG approach for CRTs proposed by 
Stephens et al. (2012). To be CAN, the DR estimator requires that either the OM or PS 
model be correctly specified regardless of the choice of the working correlation matrix. 
Interfering covariates can be ignored if either the OM or the PS is correctly specified. In 
presence of treatment-covariate interactions, if the PS is not correctly specified, covariates 
that interact with treatment on the outcome must be included in the OM. We accommodate 
these treatment-covariate interactions by modeling the OM separately for each treatment 
group. Covariates for the OM and the PS may be selected using automatic variable selection 
procedures such as a stepwise procedure, and may be at the cluster level or individual level. 

We recommend using V^W^Xi, Ai,rj w ) to ensure consistency of the IPW and the DR 
for CRTs, rather than the conventional implementation, W^ 2 {fly^V^W^ 2 (rj w ), available 
in several software packages of the weighted GEE. See Tchetgen Tchetgen et al. (2012) 
for a similar result for longitudinal data with observation-specific weights. If a working 
independence correlation structure is used, then the two implementations lead to the same 
result. When 'W 1 J 2 {r] w )V~ l W l J 2 (r] w ) and an arbitrary correlation structure is used in 
the DR, estimation of marginal treatment effect is consistent only if the OM is correctly 
specified. We provide an R package called CRTgeeDR that implements the proposed DR 
estimator. The application of our methods to data from the SAM study showed an effect of 
HIV/STI intervention on the percentage of protected intercourse (Jemmott III et ah, 2014) 
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that reached a 0.05 level of significance. Moreover, results of the analysis that distinguishes 
among different types of partners and of sexual behavior may be useful in targeting future 
interventions. Our approach allows a situation that we denoted covariate interference in 
CRTs, and thus extends ideas of adjustment of time-varying covariates in longitudinal 
responses (Pepe and Anderson, 1994; Tchetgen Tchetgen et ah, 2012). Since treatment is 
randomized at a cluster level and we consider a marginal mean model which only includes 
treatment, the covariate interference have a different implication for analysis than exposure 
interference in causal framework (Liu and Hudgens, 2014) or confounding by cluster in 
observational studies (Berlin et ah, 1999; Huang and Leroux, 2011). However, when there 
are interactions between X 7 - and Ai exposure and covariate interference are related; in this 
case, individual ij may be seen as receiving pseudo-treatment A{X j). For such a setting, our 
work may be seen as extending the notion of exposure interference in RTs to CRTs and 
is related to the work of Ogburn and VanderWeele (2014). In any case, modeling covariate 
interference may lead to substantial gains of efficiency if they predict the outcome. Therefore, 
it may be profitable to develop methods that make use of contact network information to 
inform the selection of interfering covariates. Finally, an IPW sensitivity analysis to address 
outcome MNAR as in Rotnitzky et ah (1998); Vansteelandt et ah (2007) would be useful to 
developed. 

7. Web-Supplementary Materials 

Web Appendices, Tables, Figures, simulated data and, R sources implementing the estimators 
referenced in Sections 3.1, 3.3 and 5.2 are available with this paper at the Biometrics website 
on Wiley Online Library. 
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Table 1: Descriptive statistics of outcomes, sociodemographic individual covariates and 
interfering covariates by intervention group in SAM study. 


Descriptive Statistics of the outcomes 

HIV/STI 

Control 

group 

Mean [1QRJ 

% missing 

Mean [1QRJ 

% missing 

Primary outcome for percentage of protection (Y) 

Overall 64% [26; 100J 

20.8% 

60% [22; 100] 

17.5% 

Secondary outcomes for percentage of protection (Y l ,Y' Z , V 3 and F 4 ) 

Main partner vaginal sex 61% [22; 100J 

Casual partners vaginal sex 68% [33; 100] 
Main partner anal sex 37% [ 0; 68] 

Casual partners anal sex 35% [ 0; 100] 

10.2% 

19.7% 

11.2% 

15.1% 

56% [ 0; 100] 
68% [33; 100] 
52% [ 0; 100] 
31% [ 0; 100] 

9.3% 

17.1% 

8.6% 

12.8% 



Descriptive 

Statistics of the covariates 





HIV/STI 

Control group 

p-value for association with 

Y* P(Y observed)** 

Mean [IQR] 

Mean [IQR] 

V2^0 

V3^0 

rtf ^ 0 

rtf ^ 0 

Individual covariates X yj 

Age 

26 [21; 30) 

26.5 [21; 31] 

0.41 

0.13 

0.03 

0.18 

Employment Yes 

23% 

26% 

0.04 

0.17 

0.01 

<0.001 

Married Yes 

23% 

24% 

0.05 

0.76 

0.68 

0.50 

Education Yes 

46% 

42% 

0.58 

<0.001 

0.76 

0.05 

Number of children 

1.5 [0; 2] 

1.7 [0 ;2] 

0.21 

0.12 

0.25 

0.31 

Wealth 

5.3 [4; 7] 

5.3 [4; 7] 

0.77 

0.96 

0.25 

0.54 

Social desirability 

3.4 [3.2; 3.4] 

3.4 [3.2; 3.4] 

0.87 

0.33 

0.04 

0.34 

Religiosity 

0.01 [-0.7 ;0.7] 

0.00[-0.7 ;0.6] 

0.46 

0.25 

0.07 

0.69 

HIV/STI Knowledge 

14.3 [12; 17] 

14.1 [12; 17] 

0.13 

0.93 

0.37 

0.03 

Condom Behaviors 

3.7 [3.3 ;4] 

3.7 [3.3 ;4.1] 

<0.001 

0.36 

0.16 

0.33 

Condom Knowledge 

3.1 [3; 4] 

3.1 [3; 4] 

0.41 

0.57 

0.21 

0.06 

Condom Efficacy 

3.9 [3.7 ;4.2] 

3.9 [3.7 ;4.2] 

0.01 

0.31 

0.97 

0.42 

Condom Peer norm 

3.7 [3.4 ;4.1] 

3.7 [3.4 ;4] 

<0.001 

0.71 

0.49 

0.32 

Never had HIV test 

20% 

21% 

0.61 

0.80 

0.74 

0.34 

Sexual Activity Yes 

84% 

84% 

0.71 

0.06 

0.53 

0.77 

Eating attitude 

4.2 [4 ;5] 

4.2 [3.7 ;5] 

0.76 

0.01 

0.74 

0.53 

Exercise Yes 

43% 

42% 

0.99 

0.04 

0.12 

0.46 

CAGE >= 2 

62% 

58% 

0.22 

0.41 

0.18 

0.08 

Health Knowledge 

10.8 [9; 12] 

10.6 [9; 13] 

0.51 

0.38 

0.59 

0.83 

Interfering covariates Xi. 

n,- S) = l,...,n,- 





Mean Age 

26 [25 ;27J 

27 [26 ;28J 

0.39 

0.96 

0.05 

0.10 

Mean Education Yes 

27% 

8% 

0.58 

0.61 

0.72 

1.00 

Mean Number of children 

1.6 [1.2; 2.1] 

1.7 [1.1 ;2.1] 

0.81 

0.67 

0.14 

0.59 

Mean Wealth 

5.4 [4.4 ;6.2] 

5.2 [4.4 ;6.1] 

0.45 

0.38 

0.23 

0.92 

Mean Social desirability 

3.4 [3.3 ;3.4] 

3.4 [3.3 ;3.4] 

0.16 

0.44 

0.60 

0.85 

Mean Religiosity 

0.00 [-0.1 ;0.1] 

0.00 [-0.1 ;0.1] 

0.84 

0.70 

0.18 

0.94 

Mean HIV/STD Knowledge 

14.2 [14; 15] 

13.9 [13 ; 14] 

0.37 

0.23 

0.01 

0.45 

Mean Condom Behaviors 

3.7 [3.6 ;3.8] 

3.7 [3.7 ;3.8] 

0.37 

0.40 

0.02 

0.95 

Mean Condom Knowledge 

3.1 [2.9 ;3.3] 

3.1 [2.9 ;3.2] 

0.52 

0.21 

0.15 

0.32 

Mean Condom Efficacy 

3.9 [3.7 ;4.0] 

3.9 [3.8 ;4.0] 

0.23 

0.38 

0.21 

0.58 

Mean Condom peer norm 

3.7 [3.6 ;3.8] 

3.7 [3.6 ;3.7] 

0.23 

0.52 

<0.001 

0.01 

Mean Eating attitude 

4.2 [4.1;4.3] 

4.2 [4.0 ;4.3] 

0.71 

0.15 

0.25 

0.07 

Mean Exercise Yes 

76% 

82% 

0.43 

0.53 

0.10 

0.82 

Mean CAGE>=2 

63% 

37% 

0.99 

0.79 

0.71 

0.41 

Mean Health Knowledge 

10.7 [10.5 ; 11] 

10.6 [10.3 ; 10.8] 

0.10 

0.10 

0.15 

0.73 


* Wald test for and r?!/ in the regression Y = rjg + r\\ A + rf^X + r/g AX 

** Wald test for rtf and rtf in the regression logit[P(R = 1)] = r/g 1 + rtf A + r/2 1 X + AX 
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Table 2: Analysis of effect of STI/HIV intervention on overall percentage of protected 
intercourses during the last 3 months one year after intervention (primary outcome) and 
stratified by intercourse types (secondary outcomes) in SAM study with the GEE, the IPW, 
the AUG and the DR. 


Independence (-1) Exchangeable (-E) 

{3 a SE p-value $a SE p-value 


Overall percentage of protected intercourse ( Y ) 


GEE 

IPW 

AUG 

DR 

3.751 

3.445 

5.414 

7.341 

2.419 

2.558 

1.665 

2.923 

0.121 

0.178 

0.001 

0.012 

3.738 

3.429 

5.478 

7.386 

2.361 

2.488 

1.633 

2.885 

0.113 

0.168 

0.001 

0.010 

Percentage of protected vaginal intercourse with marital partner (U 1 ) 

GEE 

5.805 

2.689 

0.031 

5.761 

2.67 

0.031 

IPW 

5.660 

2.720 

0.037 

5.626 

2.698 

0.037 

AUG 

6.550 

1.811 

<0.001 

6.518 

1.794 

<0.001 

DR 

7.254 

2.542 

0.004 

7.273 

2.50 

0.004 

Percentage of protected vaginal intercourse with casual partner (Y 2 ) 

GEE 

-0.621 

4.180 

0.882 

-0.497 

4.164 

0.905 

IPW 

-1.500 

4.182 

0.720 

-1.356 

4.17 

0.745 

AUG 

-1.191 

2.638 

0.652 

-1.121 

2.624 

0.669 

DR 

-2.103 

4.077 

0.606 

-2.018 

4.058 

0.619 

Percentage of protected anal intercourse with marital partner ( Y 3 ) 

GEE 

-0.983 

1.083 

0.364 

-0.972 

1.081 

0.369 

IPW 

-0.934 

1.087 

0.390 

-0.921 

1.085 

0.396 

AUG 

-0.951 

0.684 

0.164 

-0.954 

0.684 

0.163 

DR 

-0.835 

1.005 

0.406 

-0.819 

1.003 

0.414 

Percentage of protected anal intercourse with casual partner (U 4 ) 

GEE 

0.013 

1.201 

0.991 

-0.002 

1.204 

0.998 

IPW 

-0.003 

1.181 

0.998 

-0.019 

1.184 

0.987 

AUG 

-0.467 

0.834 

0.576 

-0.476 

0.837 

0.570 

DR 

-0.963 

1.207 

0.425 

-0.971 

1.208 

0.421 
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Tabic 3: Properties for the Doubly robust estimator (DR) compared to the GEE, the 
IPW and the AUG using the data generation mechanism from Equation 5 with covariate 
interference for the outcome and missing data generation process. Misspecihed (.MISS), 
correctly specified (.TRUE) and partially specified without treatment-covariate interactions 
(.NONE) OM and PS are investigated. Statistics for 1000 replicates are the bias compared to 
Mg = 2.0, the empirical standard errors over the replicates, the mean asymptotic nuisance- 
adjusted standard error 

and the coverage with independence (-1) and exchangeable (-E) working correlation matrix. 






Standard Error (SE) 

Coverage 



Bias 

Empirical 

Robust 

95% 


m* e 

-I 

-E 

-I 

-E 

-I 

-E 

-I 

-E 

Small sample M = 10, m - 

= (10, 20,30) with probability 1/3 each, Low correlation 



GEE (no missing) 

2.0 

0.0186 

0.0171 

0.6553 

0.6598 

0.5629 

0.5682 

93.0 

92.9 

GEE 

2.0 

-1.7186 

-1.7166 

0.5717 

0.5724 

0.5074 

0.4306 

12.8 

7.2 

IPW.PS.TRUE 

2.0 

-0.1623 

-0.1689 

1.1447 

1.1473 

0.7987 

0.8161 

83.9 

84.7 

AUG.OM.TRUE 

2.0 

-1.8142 

-1.8134 

0.4530 

0.4148 

0.8751 

0.8699 

39.4 

38.0 

DR. OM.MISS.PS.TRUE 

2.0 

-0.0127 

-0.0366 

2.7327 

2.6793 

1.4029 

1.3985 

92.0 

92.0 

DR. OM.TRUE.PS.MISS 

2.0 

0.0011 

0.0001 

0.1544 

0.1545 

0.1287 

0.1330 

86.0 

87.5 

DR. OM.TRUE.PS.TRUE 

2.0 

-0.0017 

-0.0022 

0.1881 

0.1838 

0.1413 

0.1447 

86.9 

87.4 

DR. OM.TRUE.PS.NONE 

2.0 

0.0006 

-0.0003 

0.1612 

0.1608 

0.1330 

0.1368 

85.8 

87.8 

Large sample M = 100, Hi 

= (90,100, 110) with probability 1/3 each, Low correlation 



GEE (no missing) 

2.0 

0.0042 

0.0043 

0.1156 

0.1157 

0.1155 

0.1156 

94.3 

94.5 

GEE 

2.0 

-1.7335 

-1.7321 

0.1015 

0.1013 

0.0994 

0.0994 

0.0 

0.0 

IPW.TRUE 

2.0 

-0.0113 

-0.0108 

0.2626 

0.2621 

0.2507 

0.2510 

93.5 

93.9 

AUG. TRUE 

2.0 

-1.8021 

-1.8024 

0.0694 

0.0664 

0.2556 

0.2550 

0.0 

0.0 

OM.MISS.PS.TRUE 

2.0 

-0.0089 

-0.0079 

0.3127 

0.3105 

0.3937 

0.3940 

99.3 

99.1 

OM.TRUE.PS.MISS 

2.0 

0.0013 

0.0014 

0.0259 

0.0259 

0.0256 

0.0257 

95.2 

95.7 

OM.TRUE.PS.TRUE 

2.0 

0.0013 

0.0014 

0.0284 

0.0284 

0.0285 

0.0285 

95.8 

96.0 

OM.TRUE.PS.NONE 

2.0 

0.0014 

0.0014 

0.0266 

0.0266 

0.0263 

0.0263 

95.2 

95.1 


Marginal model for the GEE: 

Ai ) — /3o + 0AAi 


OM is fitted for each treatment group Ai = a: 

OM.TRUE B(Xi, Ai = a) = 7 g + + 7 2 a XU. 

OM.MISS B(Xi, Ai = a) = 'yS+ liX2 i:i 

PS is fitted for the whole dataset: 

PS.TRUE mj (Xi, Ai) = expit (7^ + 7 % At + 71 V Xl tj + 7 f XU. + j^AiXUj) 

PS.MISS mj ( X i , Ai) = expit ( 7o M + 7 % Ai + 7^X2^) _ 

PS.NONE 7nj(Xi,Ai) - expit ( 7l f + 7 % Ai + ^ XUj + XU.) 
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Tabic 4: Sample size effect and correlation magnitude effects for data generation mechanism 
given in Equation 5 with / 3° = (1,1,1,1,1) and j3 M = (—3,1/2,1/2,1/2,1/2). Statistics 
for 1000 replicates are the bias compared to M|., the empirical standard errors over the 
replicates, the mean asymptotic nuisance-adjusted standard errors and the coverage for the 
GEE, the IPW, the AUG and the DR with independence (-1) and exchangeable (-E) working 
correlation matrix. 


Standard Error (SE) Coverage 

Bias Empirical Robust Fay’s Robust Fay’s 



Me 

-I 

-E 

-I 

-E 

-I 

-E 

-I 

-E 

-I 

-E 

-I 

-E 

Small sample M = 

10, rii = (10,20,30) with probability 1/3 each, Low correlation 






GEE 

2.0 

-1.7473 

-1.7479 

0.4351 

0.4360 

0.3963 

0.3256 

0.4559 

0.4603 

0.8 

2.3 

3.9 

4.9 

IPW 

2.0 

-1.0130 

-1.0130 

0.6793 

0.6842 

0.5538 

0.5591 

0.6735 

0.6766 

49.0 

49.2 

59.8 

59.9 

AUG 

2.0 

-1.8099 

-1.8111 

0.3371 

0.3269 

0.8362 

0.8353 

0.8834 

0.8817 

29.7 

29.1 

40.1 

39.2 

DR 

2.0 

0.0008 

0.0006 

0.1552 

0.1586 

0.1127 

0.1140 

0.1190 

0.1201 

84.8 

83.8 

86.0 

86.2 

Large sample M = 

= 100, rii = 

= (90,100, 110) with probability 

1/3 each, Low correlation 





GEE 

2.0 

-1.7335 

-1.7321 

0.1015 

0.1013 

0.0985 

0.0727 

0.0994 

0.0994 

0.0 

0.0 

0.0 

0.0 

IPW 

2.0 

-0.9955 

-0.9952 

0.1514 

0.1517 

0.1559 

0.1563 

0.1588 

0.1592 

0.2 

0.2 

0.2 

0.2 

AUG 

2.0 

-1.8019 

-1.8022 

0.0695 

0.0664 

0.2556 

0.2550 

0.2569 

0.2563 

0.0 

0.0 

0.0 

0.0 

DR 

2.0 

0.0016 

0.0017 

0.0265 

0.0265 

0.0262 

0.0263 

0.0264 

0.0264 

95.1 

95.0 

95.1 

95.2 

Small sample M = 

10, rii = (10,20,30) with probability 1/3 each, High correlation 






GEE 

2.0 

-0.0086 

-0.0086 

0.5265 

0.5314 

0.4701 

0.4721 

0.5651 

0.5657 

88.5 

88.4 

92.9 

92.7 

IPW 

2.0 

-1.0221 

-1.0229 

0.7026 

0.7083 

0.5776 

0.5829 

0.7015 

0.7044 

52.4 

52.2 

62.2 

61.5 

AUG 

2.0 

-1.7985 

-1.7987 

0.5058 

0.5084 

0.8748 

0.8727 

0.9243 

0.9209 

35.8 

35.8 

45.1 

45.5 

DR 

2.0 

0.0098 

0.0062 

0.4328 

0.4407 

0.2469 

0.2480 

0.2607 

0.2614 

77.4 

77.7 

79.7 

79.6 

Large sample M = 

100 , rii = 

(90, 100, 110) with probability 1/3 each, High 

correlation 





GEE 

2.0 

-1.7325 

-1.7312 

0.1145 

0.1141 

0.1121 

0.0753 

0.1132 

0.1132 

0.0 

0.0 

0.0 

0.0 

IPW 

2.0 

-0.9945 

-0.9940 

0.1618 

0.1620 

0.1652 

0.1656 

0.1682 

0.1686 

0.2 

0.2 

0.2 

0.2 

AUG 

2.0 

-1.8014 

-1.8017 

0.0787 

0.0761 

0.2587 

0.2581 

0.2600 

0.2594 

0.0 

0.0 

0.0 

0.0 

DR 

2.0 

0.0029 

0.0032 

0.0609 

0.0610 

0.0590 

0.0590 

0.0593 

0.0593 

94.7 

94.6 

94.7 

94.6 


Marginal model for the GEE: 

/l(/h A) — A) + /UA 

OM in AUG and DR is fitted for each treatment group A = a using a stepwise regression: 

B(Xi,Ai = a) = stepwise(Xlij, X2ij, X3ij, Xli., X2i., X3i.) 

PS in DR and IPW is fitted for the whole dataset using a stepwise regression: 

loQit(^7V{j (.A i , Ai ) ) — stepwise(Aj , X l,.y . X 2^j , X 3,j , X U. , X 2 j. , X 3 1 , ) 









Marginal treatment effect in CRTs with missing outcomes 
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Tabic 5: Simulation of the scenario described in Equation 6 mimicking the SAM study data. 
Statistics for 1000 replicates are the bias compared to M the empirical standard errors 
over replicates, the mean asymptotic nuisance-adjusted standard error, and the coverage for 
the GEE, the IPW, the AUG and the DR with independence (-1) and exchangeable (-E) 
working correlation matrix. 


Standard Error (SE) Coverage 

Bias Empirical Robust Fay’s Robust Fay’s 



Ml 

-I 

-E 

-I 

-E 

-I 

-E 

-I 

-E 

-I 

-E 

-I 

-E 

Small sample M = 10, 
GEE 5.73 2.214 

m = (10, 20, 30)with probability 1/3 each 
2.213 1.330 1.329 1.829 1.848 

1.359 

1.363 

89.4 

88.7 

59.4 

59.7 

IPW 

5.73 

0.536 

0.537 

1.333 

1.333 

1.214 

1.214 

1.470 

1.471 

86.5 

86.4 

92.1 

92.1 

AUG 

5.73 

0.173 

0.173 

0.973 

0.973 

0.878 

0.878 

0.925 

0.925 

88.6 

88.6 

89.9 

89.9 

DR 

5.73 

-0.104 

-0.104 

1.102 

1.101 

0.932 

0.931 

0.982 

0.982 

90.3 

90.3 

92.0 

92.0 

SAM-like 

GEE 

sample M = 
5.73 2.347 

50, m = 
2.343 

(20, 30, 30) with 
0.308 0.308 

probability 1/3 each 
0.532 0.466 0.308 

0.309 

0.0 

0.0 

0.0 

0.0 

IPW 

5.73 

0.622 

0.623 

0.303 

0.303 

0.317 

0.317 

0.323 

0.323 

50.7 

50.7 

52.1 

52.1 

AUG 

5.73 

0.215 

0.215 

0.222 

0.222 

0.230 

0.230 

0.232 

0.232 

85.1 

85.1 

85.2 

85.2 

DR 

5.73 

0.037 

0.026 

0.259 

0.260 

0.252 

0.253 

0.254 

0.255 

94.6 

95.3 

94.6 

95.4 


Marginal model for the GEE: 

M(/h -A) — ho + /A A 

OM in AUG and DR is fitted for each treatment group using a stepwise regression: 

B(Xi,Ai = a) = stepwise(EMP ii ,MARI ij -, AGE ij -,REL ij -,CAGE ij - ! HIV ii ,CDM ij -,Xl ii ,X2 i:# ,X3y) 

PS in IPW and DR is fitted for the whole dataset using a stepwise regression: 

logit(TVij(Xi, Ai)) = stepwise(Ai, EMPij, MARU,-, AGEij, RELi-,-, CAGEij, HIVij, CDMij, Xl, 3 -, X2ij,X3y) 







